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For a QCD chiral critical point to exist, the parameter region of small quark masses for which the 
finite temperature transition is first-order must expand when the chemical potential is turned on. 
This can be tested by a Taylor expansion of the critical surface (m„ ^/,mj)c(/x). We present a new 
method to perform this Taylor expansion numerically, which we first test on an effective model of 
QCD with static, dense quarks. We then present the results for QCD with 3 degenerate flavors. For 
a lattice with Nt —4 time-slices, the first-order region shrinks as the chemical potential is turned 
on. This implies that, for physical quark masses, the analytic crossover which occurs at jU = 
between the hadronic and the plasma regimes remains crossover in the -region where a Taylor 
expansion is reliable, i.e. pL<T . We present preliminary results from finer lattices indicating that 
this situation persists, as does the discrepancy between the curvature of Tc{mc{pL —Q),li) and the 
experimentally observed freeze-out curve. 
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Figure 1: Schematic phase transition behaviour of Nf = 2+ \ flavor QCD for different choices of quark 
masses (ot„ ;/,OTs) at /i = [|i]| {left), and numerical determination of the chkal critical line [Q] (right). 



1. Introduction 



The experimental determination of the QCD phase diagram is underway via extensive heavy- 
ion collision programs. At the same time, much effort is being devoted to its theoretical determi- 
nation via numerical lattice simulations. In the latter case, one can ask the more general theoretical 
question: what is the phase diagram of 2+ 1 -flavor QCD in the (m„,i,ni.s,/i„ j,/Xi,r) parameter 
space? 

For iJLu4 = M.V = 0, theoretical expectations are summarized in Fig. 1 (left). In the chiral (m„ ^/ = 
fUs = 0) and the quenched {niuj = = °°) corners, order parameters probing the breaking of the 
chiral and the center symmetries exist, and the symmetry breaking or restoring transitions are first- 
order. For intermediate quark masses, simulations show a crossover. This leads to the existence of 
two lines of critical points, both in the 2>d Ising universality class, delimiting the first-order regions. 
In the quenched corner, the critical line has been studied in [|2|, ^. In the chiral corner. Fig. 1 (right) 
shows the result of [Q]. Good agreement with expectations was found, including consistency with 
a tricritical point (m„^^ = Oj/ti^ = m'™) with a rather heavy mass m^"'^ ~ 2.87^. The physical point, 
marked by an X, lies in the crossover region. These results were obtained with standard staggered 
fermions on an A/^r = 4 (a ~ 0.3 fm) lattice. One task is now to quantify cutoff effects and extrapolate 
to the continuum limit. This extrapolation has been performed for the physical point [§], confirming 
that it remains in the crossover region. Our preliminary results on an A', = 6 lattice, consistent with 
those of ||^ and with earlier indications from improved actions on A'^r = 4 |^, are presented Sec. 3^. 
Cutoff effects are large, and the transition becomes much weaker on a finer lattice. 

The next issue is to determine the effect of a baryonic chemical potential. For physical quark 
masses, it is expected that the transition becomes stronger with }JL, so that it turns from a crossover at 
/I = into a first-order transition, thus defining the QCD critical point T^) where the transition 
is second-order. The fermion determinant becomes complex when /i 7^ 0, and the ensuing "sign 
problem" forbids standard Monte Carlo sampling. Nevertheless, a determination of this critical 
point has been performed, using the same staggered action and Nt = 4 lattice spacing as us [||]. 
While the cutoff error is likely to be comparable to that at /i = 0, the numerical method used. 
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Figure 2: For physical m^) quark masses, the finite-temperature "transition" at /i = is really a 

crossover. As pL is turned on, this crossover becomes a genuine phase transition at the QCD chiral critical 
point (left), provided the region of first-order transitions at small quark masses expands with jU. If not {right), 
there is no QCD chiral critical point. The curvature of the critical surface dmc/djj.^ at jj. —0 distinguishes 
between these two possibilities for small jj.. Note that, for heavy quarks, the first-order region shrinks when 
jj, is turned on 

namely reweighting, introduces new possible systematic errors. Contrary to the well-known multi- 
histogram reweighting where one interpolates between Monte Carlo results obtained at several 
values of the coupling, here one extrapolates results obtained at a single value of the coupling, n = 
0, to /I 7^ 0. While this procedure is in principle exact for infinite statistics, the practical question 
is whether the finite Monte Carlo sample contains the relevant information at the extrapolated 
coupling. This difficulty is known as the overlap problem. Thus, a crosscheck of the results of 
using a different approach seems worthwhile. 

To avoid all difficulties caused by a potential lack of overlap, we set a more modest goal, and 
study the effect of an infinitesimal chemical potential. More precisely, we consider the critical 
surface, swept by the chiral critical line of Fig. 1 (right) as a function of /i, in terms of a Taylor 
expansion in {ji/TY about jJL =0, (for the Nf = 3 case to keep the notation simple) 



The sign of the first Taylor coefficient c\ is of crucial relevance to the QCD critical point. Since 
for /I = the physical point is in the crossover region, the first-order region must expand with 
jX (Fig. 2 left). If instead the first-order region shrinks (Fig. 2 right), then the crossover persists 
and there is no QCD critical point unless another critical surface, topologically unrelated to the 
chiral critical surface we consider, is present in the phase diagram. The simplest corresponding 
(/i,r) phase diagram in the two cases is depicted in Fig. 3, for a quark mass below and above 
mc(0). In the "exotic" scenario (right), the first-order line present for small quark mass turns into 
a crossover as the mass increases, starting with large pt. If m > mc(0) as in the physical situation, 
the "transition" is simply a crossover, for all /I's, until different physics take over. Distinguishing 
between these two scenarios requires full knowledge of the critical surface. However, by measuring 
the first Taylor coefficient, we can determine the behaviour of the critical surface for jU/r < 1, i.e. 
in the region pts ^ 500 MeV where experimental searches are considered. 




(1.1) 
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Figure 3: Simplest phase diagram in the (/i , T) plane, for the two cases of Fig. ^(/e/f and right), as a function 
of the quark mass {middle). Red and blue lines indicate first-order transition and crossover, respectively. 



In our first study we performed simulations at imaginary /i and determined mc{}JL = ipLi), 
then fitted our results by a truncated Taylor expansion which can be trivially continued to real 
/I. While this method appears to work very well for the pseudo-critical line T^{pL) [ p^ , 11], we 



faced two related difficulties when computing mc{)x): the signal is very noisy, and this makes the 
determination of the systematic error due to truncation of the Taylor series difficult. For the case of 
three degenerate flavors, = m„ d, /i^ = /Xj, ^, we performed quadratic and quartic fits in n/T. The 
quartic term was statistically insignificant, so we set it to zero and obtained for the leading term 

which favors the "exotic" scenario Fig. 2 (right), but not even at the 2a level. Note that including 
the quartic term in the fit (Table 2, line 2 of Ref. [||]) changes the leading coefficient to -2.6(1.2). 
The lack of compelling numerical evidence motivated us to improve our statistics and our numerical 
methods. 

Here, we compare three methods to obtain the derivative ^|;U=o- 

A. Analytic continuation from imaginary /i 

B. Direct measurement of the derivative at /i = 

C. Noisy reweighting to small imaginary /i 

In all three methods, the same criterion is used to determine the pseudo-critical couplings in a 
finite volume: the Binder cumulant B4 = -^^^Ip-, with 5X = X — {X) and X = yy, takes the 
value 1.604 characteristic of the 3d Ising universality class at the critical point. In a finite volume, 
B4 is an analytic function of (m,r, /x), and T is fixed to the pseudo-critical temperature Tc{m,}jL) 
by requiring, for instance, that the susceptibility be maximum. We use the equivalent 
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prescription {{5X)^) = 0. The two constraints {S4 = 1.604, = 0} define a line in the 

{T,m,ix) space. We want to extract the curvature ^ of this line at = 0. 

Method A is the one used in The results of different simulations at various imaginary 
values of pL are fitted with a truncated Taylor expansion about /i = 0. One difficulty comes from 
the systematic error associated with the order of the truncation. Another comes from the final 
statistical error, which in [Q] amounted to about half the signal. 

Method B is reminiscent of the Bielefeld-Swansea approach [11], where the derivatives of 
the free energy with respect to jx/T are expressed as expectation values of non-local operators, 
which involve traces of inverse powers of the Dirac operator. In our case, the derivatives of B4 
and {{5\j/Yy) with respect to m,T and jj. can be expressed as expectation values of complicated 
operators, which can be measured in a single simulation at = 0. 

Method C, to our knowledge, is a new application of reweighting, where the reweighting 
factors are not evaluated exactly but estimated by a stochastic method. This noise does not prevent 
reweighting provided it is unbiased. We can apply this strategy to the results of /X = simulations 
and reweight them to imaginary /i, monitoring the change in B4 and {(dxjfY)^)- Keeping this 
imaginary /i very small has two advantages: (/) the reweighting factor remains close to 1 and the 
overlap of the jJ. =0 Monte Carlo ensemble with the target At 7^ ensemble is guaranteed; (//) 
fluctuations in the two ensembles are strongly correlated, and cancel in the observables. 

In Section 2, we first test methods B and C on a Potts model representative of the heavy-quark 
limit of QCD, for which we have previously determined the critical line for both imaginary and real 
/I [^. All methods agree, with a clear advantage to method C for its simplicity. Then in Section 
3, we compare methods A and C for Nf = 3 QCD, on a coarse 8^ x 4 lattice. Again, methods A 
and C agree, leaving no doubt that for this system the first-order region shrinks when the chemical 
potential is turned on. Finally, in Section 4 we present some preliminary results on a finer 18^x6 
lattice, and discuss the implications and limitations of our findings. 



2. Potts model 

2.1 An effective description of dense static quarks 

When the quark mass m becomes infinite, quarks become static and a quark source is a 
Polyakov loop. The canonical partition function of QCD with n static quarks and fi static anti- 
quarks is [ ]T2| ] 

/■ 1 1 - m 

Zn,n= ^A-^[A]"-^*[A]"exp{-in + fi)-)exp{-Sg[A]) , (2.1) 

where ^[A\ is the Polyakov loop integrated over space and Sg the gauge action. The grand- 
canonical partition function is then simply 

Z(At) = £^'("-")^Z„,,= /^Aexp(-5,[A]+e-^<I>[A]+^'-'?^<I>*[A]) . (2.2) 

n,n 

Simplifying the gauge action to a Potts interaction between neighbouring Polyakov loops, one 
obtains the partition function of a 3(i ^ = 3 Potts model in an external field: 

Z{K,m,fi)= ^ exp[-K'£5^^^),^^jf+i)+Y,ih^{x)+h'^*{x))] , (2.3) 

mx)} Lx X 
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Figure 4: Phase diagram of the Potts model, as a function of m/T and {jj./T)^. For heavy masses, the 
transition is first order. It turns into a crossover for light masses. The critical line can be determined for 
imaginary jj, (red points), showing the expected Z3 transition dXfx/T^ i^, or directly for real /i, because the 
sign problem is mild. The line is analytic at /i = 0, but shows some curvature. Methods B and C measure its 
slope at /X = 0. 



where is the Potts (actually Z3) spin exp{ik{x)^), h = hme^^ ,h' = h^e^^ ,hm = e^"\jl = j 
and m = J. Contrary to [ p^ who consider the limit m, /I — > 00, — ^) finite, and drop the last term 
corresponding to antiquarks, we keep the complete expression in order to study the {rh,p.) phase 
diagram. Because h' / h* when /i 7^ 0, the action becomes complex and a sign problem appears. 
However, it is very mild, so that simulating with the real part of the action and reweighting with the 
imaginary part works efficiently and reliably. Two phases are present: confined/disordered at small 
K, deconfined/ordered at large K. The phase transition is first-order for the ordinary Potts model, 
i.e., in the absence of any external field h = h' = 0. This case corresponds to m — > 00. 

Decreasing m turns on the magnetic field. This weakens the phase transition, so that it becomes 
a crossover for sufficiently large h^, i.e., sufficiently small fh = m/T . The critical value of m corre- 
sponding to the end of the first-order region depends on the chemical potential /i. The critical line 
tnc{jJ.) has been determined in [Q], for both real and imaginary chemical potentials. Fig. ^ sum- 
marizes the results of [^. It represents the qualitative behaviour of QCD with heavy quarks near 
the critical line in the upper right corner of Fig. 1 (left), because the symmetries and the infrared 
degrees of freedom are the same. Given Tc ~ 270 MeV, the value of mc{0) even lies in the 1-2 GeV 
range estimated from full QCD simulations |Q, |3|]. Fig. 4 shows that, as the chemical potential is 
turned on, the first order region shrinks as indicated in the quenched corner of Fig. 2 (right). 

As a function of p.^, the critical line is analytic at /I = 0, so that the real-/i dependence can be 
reconstructed by analytic continuation of imaginary-/! simulation results. However, some curvature 
is visible, particularly as the imaginary p. approaches the Z3 transition point fl = ij. A global fit of 
all our data requires a fourth-order polynomial in jl^: 

^ = 8.273 + 0.585 (^)'-0.174(^)Vo.l6o(^)'-0.07l(^)' , (2.4) 

where the first coefficient is 0.585(3) including its fitting error. This is essentially the result of 
method A (although we have actually augmented the imaginary-/! data with real-/! data). One 
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may be concerned that a similar curvature occurs in QCD as well, which would make it difficult to 
analytically continue data from a few discrete values of imaginary pL. Therefore, the object of this 
study of the Potts model is to see how well the first coefficient of the polynomial, 0.585(3), can be 
reproduced by methods B and C which involve simulations at /X = only. 

2.2 Method B 

The order parameter probing the Z3 symmetry is the magnetization M = Y,x^{^)- Therefore, 
we determine the critical line by requiring 



^4-^^^^ = 1.604, {{5M)') = , (2.5) 



with 5M = M— (M). Both expressions are functions of (K:,m,/I), so that the two equations deter- 
mine the critical line rhcifl) in this 3-parameter space. The two equations can be Taylor expanded 
about the point {Kcifi = 0),mc{jX = 0),0) of this line, yielding at lowest order 

dB4=Adm + Bdll^ + CdK, d{{5Mf)=A'dm + B'dll^ + C'dK . (2.6) 

Staying on the critical line implies that the changes in B4 and in ((5M)^) are both zero. From the 
resulting two linear equations, one can eliminate dK, and obtain 

dm _ BC'-B'C 

dp^~~AC'-A'C ' ^^^^ 

which is the desired slope of the critical line at /i = 0. The coefficients A,B,C,A' ,B' ,C' are the 
following expectation values 

A = ^\f,=o = -2hU5MY^{5M') 



m 

B = -^\f,=o = ^{5MY^[2hnr{dM')+hl{{dM'-{dM')){^-^*'^' 

-hi {5mY^ {5M^) {{5M^ - {5M^)) 
C = ^|/i=o = {5mYHSM^SE)-2{5mY^5M^){5M^5E) 

A' = ^|,.o = -2/..(5M^) 
dm 

B' = \ ^^f^2^ Im-o = \[2h^{m^-H8M^)8M)M)+hl{{8M' - ^{8M^)8M){^ -^*f)] 
C = ^^|^|/i=o = {8M^8E)-3{8M^){8M8E) (2.8) 

(j IC 

where 8E = E — {E). Thus, all 6 coefficients can be measured during a Monte Carlo simulation at 
the jJ. = critical point, i.e., where K and m have been tuned to satisfy eq.(p^). We have performed 
such measurements on a 72^ lattice, of the same size used to obtain Fig. 4. Substituting into 
eq.(2.7), and performing a jackknife bin analysis, we obtain drh/dpL^ = 0.593(8), which agrees 
very well with the earlier result from method A. Note again that method B is insensitive to the 
curvature of the critical line, unlike method A. 



1 d'^BA , 1 
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Figure 5: Effect of a small change in the magnitude {left) and the orientation {right) of the magnetic field 
h — h,n exp{ijli) on the Binder cumulant B4{Re{M)) in the Potts system. 



2.3 Method C 

The previous method enforced algebraically that the Binder cumulant B4 stay constant under 
an infinitesimal change dm, dpL^. Instead, one can measure the change in B4, AB4, under a small 
variation Am or A/i^, thus estimating the finite differences AB4/Am and AB4/A/i^. For sufficiently 
small variations, these discrete differences will approach the derivatives 884/ dm and 984/ dji^. 
Finally, once the gradient of B4 is known, the direction which keeps 64 constant is given by 

Am dB4 dB4 
A/i-^oA/l^ dpL^ dm 

The second equation in ( ^^ ) is not needed here. The measurement of B4 at the shifted couplings 
{Mc + Am, n = 0) or {nicAjl) is performed by the usual analysis, which automatically tunes fc to 
satisfy {{5M)^) = 0. 

To measure B4 at the shifted couplings, one could perform new Monte Carlo simulations. But 
this is not necessary. Because the shift in the couplings is very small, it is adequate and safe to 
use the original Monte Carlo ensemble and simply reweight the results in the standard way [jPI]. 
Moreover, to avoid complex weights one reweights to imaginary p., which simply introduces a sign 



flip in eq.(2.9). The reweighting factors remain real positive, and close to 1. 

With reweighting, the fluctuations in the original and the reweighted ensembles are strongly 
correlated. This can be turned into a virtue, as these correlated fluctuations drop out of our observ- 
able, which is the change in B4, rather than B4 itself. 



8 



QCD critical point 



Philippe de Forcrand, Owe Philipsen 



This procedure is illustrated Fig. 5. The top row shows the change in B4 under a change 
in hfn = exp(— m) (left), and under a change in imaginary /i (right). One observes the expected 
dependence, respectively linear and quadratic. Note that the changes in B4 are small, ^(10-2), and 
measured very accurately - much more accurately than B4 itself. The bottom row of Fig. 5 shows 
the change in B4, divided by the change in /i,„ (left) or by (/I)^ (right). If the Taylor expansion could 
be truncated to leading order, the data would be constant. Instead, one sees the small influence of 
the next Taylor order. A fit to (constant + linear) gives the desired partial derivatives, marked by a 



black circle, which can be substituted in eq.(2.9). The resulting slope is 0.589(7), again consistent 
with the other two methods. 

Note that the statistical errors on the various points Fig. 5 are extremely correlated with each 
other, so that a jackknife bin analysis is required to obtain reliable errors on the final slope. Note 
also that there is a broad optimum for the shift in the couplings: too small a shift produces too 
small a change in 64, and the estimates of the derivatives approach 0/0; too large a shift introduces 
a systematic error from higher-order Taylor terms and a potential overlap problem. 

The errors from methods B and C are similar. This is normal, since these two methods make 
use of the same Monte Carlo data and extract the same information. So the preference for method 
C in this case comes from its simplicity: there is no need to measure the relatively complicated 



observables eqs.(2.8). This will become a more serious issue in the case of QCD. 



3. % = 3 QCD 

In [Q], for Nf = 3 QCD, with standard staggered fermions on anNt =4 lattice, we determined 
the critical quark mass at /i = 0: am^ = 0.0263(3), and proceeded to Taylor expand the pseudo- 
critical coupling j3c and the Binder cumulant of \ff\if' 

I5c{aij.,am) = ^ Cki{aii)^'^ {am — amo)' (3.1) 

k,i=o 

B^{am,a^) = 1.604 + Z7io \am — amQ — c\{a}x)^^^+b2Q{am — amQf^ 

-bm[{c'2-c\C){aixf +C{am-arrfQ){aixf]+--- , (3.2) 

from which one can extract the variation of the pseudo-critical temperature and of the critical quark 
mass with the chemical potential. 

3.1 Recall: method A 

In [Q], we performed independent simulations at different values of the quark mass and imag- 
inary chemical potential /i = //X/, which we then fitted with the Taylor expansions above. 

For the pseudo-critical temperature, a leading order fit was satisfactory, yielding for the /x^ 
dependence cio = 0.781(7). At subleading order, there was no evidence for a cross-term {am — 
amQ){aiJ.)^, and the {ajj.)'^ term was barely statistically significant. Including it in the fit yielded 
cio = 0.759(22) (see Table 1 of [@]). Both fits are shown in Fig. 6 (left), by the narrow blue and 
broader green error band, respectively. 



For the curvature of the critical surface, the coefficient c[ eq.(P^) encodes the relevant infor 



mation: to keep B4 constant, one must satisfy = c[ . Again, no cross-term C was visible. 
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and fits including ji^ only, or ji^ and ji'^ terms were performed, resulting in the narrow blue and 
broader green error bands Fig. 6 (right), corresponding to Table 2, lines 2 and 3 of [Q]. They suggest 
a negative value for c\ , but do hardly more than that. 

3.2 Method B 



The Taylor coefficients in eqs.(3.1,3.2) can be expressed as expectation values to be measured 



at jU = 0. We wrote down these expressions, analogous to eqs.(2.8). But we did not implement 
a program to measure these operators, for several reasons. First, the programming effort is non- 
trivial. For instance, the trace of the 5^'' inverse power of the Dirac operator must be evaluated to 
obtain the derivative of {5\j/Yy with respect to the quark mass. Moreover, important cancellations 
will take place among the various contributions, leading to further difficulties with optimizing the 
number of noise vectors to be used as stochastic estimators of the various traces. Finally, we 
realized, from the Potts test case, that method C makes use of the same information contained in 
the jU = Monte Carlo ensemble, and gives the same output as method B with less effort. 

3.3 Method C 

As in the Potts case, one can shift very slightly the quark mass and the chemical potential, 
and reweight the n = Monte Carlo ensembles to these shifted couplings. The effect of a shift in 
the quark mass was already measured in with sufficient accuracy, so we were interested in the 
effect of a small chemical potential, taken as imaginary to preserve positivity of the weights. The 
difference with the Potts case is that the reweighting factors, 

for each configuration {U} were not computed exactly, but only estimated as 

p([/,At2,Aii) = (expf-|^-^//^([/,At2)4^+^^/^t/,Mi)r]P + |r]|2)\ , (3.4) 



where T] is a Gaussian random vector. Since Nf = 3 in our case, the fractional powers of the Dirac 
operator were approximated to high precision by a ratio of polynomials. To reduce the variance, we 
actually formed uncorrected estimators for -y/p(?7,/i2,Mi) and multiplied them together. Further, 
p{U,lJ.2,0),lJ.2 > Hi was constructed as p{U,jX2,i^\) x P{U,IjLi,0). In this way, highly correlated 
reweighting factors were obtained for 6 values of = //i/, with ajj.] ranging from 0.01 to 0.1. 

Crucially, the noise in the reweighting factors does not prevent the standard application of 
reweighting, because these factors enter linearly in the numerator and denominator of the reweighted 
expectation value {W) = Difficulties arise only when one has to form {W)'^, because the 

same stochastic estimator for configuration / is used k times and the statistical error, when taken 
to an even power, introduces a bias. This is a relevant concern, since we want to measure the 4''' 
cumulant of yy- However, this bias is of order tff{l/N'^^^) and disappears as the Monte Carlo 
sample size grows. To test for such bias, we magnify it by multiplying every estimated weight 
by a random number drawn uniformly in ]0,2[. The results of the analysis, for the derivative cio 
of the pseudo-critical coupling eq.(0|) and for the coefficient c\ in the Binder cumulant expansion 



eq.(3.2), are presented Fig. 7. They are to be compared with Fig. 6, where the original weights 
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Figure 6: Upon introducing a small imaginary chemical potential jU/, the pseudo-critical coupling jS^ and 
the Binder cumulant £4(1//^/) vary slightly from their /i/ = values. The change /A/i^ is shown as a 
function of (a/X/)^, for (? = /3( (left) and for 0" — B4 {right). The error bands correspond to the leading-order 
(blue) and subleading-order (green) fits from Ref.[0]. 
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Figure 7: Check of possible systematic error. Same as Fig. but the reweighting factors for each configura- 
tion have been multiplied by a random number uniformly distributed in ]0 : 2[. This unbiased noise, similar 
to the noise in the reweighting factors themselves, does not bias the result. 



were used. No bias is visible in Fig. 7, leading us to conclude that no bias is present in Fig. 6 either. 
Note that our sample size is very large: we analyzed about 5 million S^* x 4 configurations at 21 
j3 -values. This large statistics was made possible by using the Grid at CERN. The actual running 
time was less than two weeks. 



3.4 Comparing methods A and C 

The final results from method C (Fig. 6) are 



CIO = 0.746(1) or fr^TT = 1 " 0-637(1) f-^V + •• • (3-5) 

ro(m^,0) \nTQj 

for the pseudo-critical temperature, using the two-loop j8-function to convert to physical units, and 



11 



QCD critical point 



Philippe de Forcrand, Owe Philipsen 























=0.00 H 






an 


=0.14 






y' an 


=0.17 






an 


=0.20 


-A 




an 


=0.24 


▼ 



0.02 



0.024 



0.028 
am 



0.032 



0.036 




m=0.020 ' — H- 

m=0.025 ♦ 

m=0.0265 A; 

m=0.0275 ;f 

m=0.030 - ♦ 

^ , , m=0,034 ; 




0.0125 0.025 0.0375 0.05 0.0625 

>2 



(am) 



Figure 8: Direct measurements of the Binder cumulant (left), updated from Ref. [|4|. Subleading depen- 
dence on m becomes visible. Comparison with the reweighting approach as a function of /X^ (right) shows 
remarkable agreement, both in the leading and subleading terms. 



for the curvature of the critical surface, where the value of biQ eq.(3.2) was taken from [Q]: biQ ^ 
13.6. A subleading term would be visible as a slope in the fits Fig. 6. Indeed, a small effect is 
visible on the right, corresponding to = —2.5 ± 1.2. 

These results are consistent with those of method A, provided subleading terms are included 
in fitting the latter (green error bands in Fig. 6), even if they are statistically almost unconstrained. 
This illustrates the difficulty of estimating the systematic error of truncating the Taylor expansion 
used in fitting. This difficulty is eliminated in our new method, which in addition is about 100 
times more efficient. 

However, method A gathers statistics over a broader range of chemical potentials, and with 
updated statistics exceeding 15 million configurations, now allows us to clearly identify subleading 
Taylor terms. They are visible from the curvature of the fits of the Binder cumulant as a function of 
the quark mass Fig. 8 (left). The complete ansatz eq.(^^) was used, now yielding c\ = —0.074(28) 
and C2 = —1.0(5) with a of 23 for 21 d.o.f. Note the consistency of methods A and C, also for 
the subleading order whose sign contributes to further shrinking the region of first-order transitions. 

Fig. 8 (right) shows the same results, after subtraction of the fitted mass dependence, as a 
function of (a/^/)^. The fit is shown by the lower parabola. Now, the results of Fig. 6 (right) 
(leading and subleading terms) are shown in the same figure as the upper parabola. The agreement 
between the two independent methods is remarkable, given that method C only probes the region 
{apLiY ^ 0.01 where agreement is near-perfect. 



3.5 Towards the continuum limit 

Of course, cutoff errors on our A^, = 4 (a ~ 0.3fm) lattice can be large, and it is essential to 
perform a continuum extrapolation. To this end, we are pursuing our project on Nt = 6 lattices, and 
present some preliminary results Fig. 9. 

The left figure illustrates cutoff effects on the critical bare quark mass, mj,, corresponding to 
a second-order transition at /i = in A/y = 3 QCD. One can see that the quark mass, expressed 
in units of the temperature, must be reduced by a factor ~ 5 on A/f = 6 lattices. A similarly large 
effect is present in the resulting pion mass, m^, measured at zero temperature for quark mass 
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Figure 9: Left: determination of the critical quark mass in the Nf ~3 theory. The bare quark mass decreases 
(in units of Tc) as the continuum limit is approached. The corresponding pion mass (measured at T = 0) also 
decreases. Right: preliminary result for the curvature of the pseudo-critical line for Nt =6. In physical units, 
the curvature is smaller than for A^, = 4. 



m^. The ratio m'^/r,. decreases from 1.680(4) (A', = 4) to 0.954(12) (Nt = 6), so that a naive 
extrapolation would give ~ 0.4 in the continuum! This very large cutoff effect is consistent with 
earlier indications and with a new study all suggesting that the transition becomes much 

weaker in the continuum limit. Note that the cutoff effect on the hadron spectrum is comparatively 
mild, so that the net effect of a finer lattice is to dramatically push the critical surface Fig. 2 toward 
the origin, while leaving the physical point untouched. Thus the gap between the critical surface 
and the physical point widens, pushing the critical point in Fig. 2 (left) to larger values of He- 

Our second preliminary result. Fig. 9 (right), shows the curvature of the pseudo-critical cou- 
pling d^c / d{a}xY\^=Q for the critical quark mass /Mq. The error band corresponds to the Nt = 4 
study. The trend is for d^c/d{aiJLY to be smaller for A/f = 6, while if we use the two-loop /3-function 
to convert to physical units, one should observe Jj8c/fi?(a/i)^ Nf. Instead of increasing by (6/4)^, 
our measured value seems to decrease. Now, for Nt = 4 already, the estimated curvature of the 
pseudo-critical line Tc{m^^,iJ.) was about 3 times less than that of the experimental freeze-out curve 
[15]. These two curves appear to become more clearly separated as a ^ 0, which also reduces Mq. 



4. Discussion 



Our findings including next-to-leading terms onN =4 predict the "exotic" scenario of Fig. 2 
(right): the region of first-order transitions shrinks as the chemical potential is turned on, so that the 
chiral critical surface does not intersect the physical line, and there is no chiral critical point in 
QCD. This statement, which goes against conventional wisdom, gets qualified by a number of 
systematic errors. Like [4-8,10,11,14], we use staggered fermions with the rooting trick, which 
is potentially unsafe for very light quark masses. Next, the curvature of the critical surface varies 
with the cutoff, and presently the rate of this change is unknown. It also changes from the Nf = 3 
case presented here to the Nf = 2 + 1 theory, although we found in [Q] that for Nt =4 the sign of 
the curvature remains unchanged. Finally, once ~ higher order terms may become relevant. 

Despite these caveats, we believe the qualitative picture to be robust. If the continuum critical 
quark mass can be Taylor expanded as per eq.(|l.l|) with coefficients ^(1), the critical surface in 
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Fig. 2 rises "almost vertically" no matter the sign of its curvature, and a critical point at small 
|/i/r| implies a fine-tuning of the physical quark masses, so as to be very close to the critical line 
at /X = 0. Such a fine-tuning seems unnatural, and indeed the = critical line seems to recede 
considerably in the continuum limit, now requiring a large curvature of the opposite sign to what 
we observe in order to accomodate a critical point at |ju/r| < 1. 

Finally, the object of our study is the chiral critical surface. Our findings do not exclude 
additional critical structure due to non-chiral physics causing the phase transitions Fig. 3, bottom 
right. 
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